Period doubling, two-color lattices, and the growth of swallowtails in Bose-Einstein 

condensates 
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The band structure of a Bose-Einstein condensate is studied for lattice traps of sinusoidal, Jacobi 
elliptic, and Kronig-Penney form. It is demonstrated that the physical properties of the system are 
independent of the choice of lattice. The Kronig-Penney potential, which admits a full exact solution 
in closed analytical form, is then used to understand in a novel way the swallowtails, or loops, that 
form in the band structure. Their appearance can be explained by adiabatically tuning a second 
lattice with half the period. Such a two-color lattice, which can be easily realized in experiments, 
has intriguing new physical properties. For instance, swallowtails appear even for weak nonlinearity, 
which is the experimental regime. We determine the stability properties of this system and relate 
them to current experiments. 
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I. INTRODUCTION 

There has been great interest in the study of Bose- 
Einstein condensates (BEC's) from both the experimen- 
tal and theoretical perspective since they were first cre- 
ated in 1995 In particular, the examination of 
BEC's in periodic potentials in the supcrfluid phase has 
yielded many intriguing phenomena, such as gap soli- 
tons 0, IU, the appearance of swallowtails or loops in 
the band structure, pulsed atom lasers and demonstra- 
tions of their phase coherence H Q, and matter- wave 
diffraction [8j|. BEC's in periodic potentials, unlike other 
solid state systems, have the advantage that the lattice 
geometry and interatomic interactions are highly control- 
lable [3, E3 ■ We examine the mean field Bloch states of a 
BEC in a one-color and two-color potential for arbitrary 
interaction and potential strengths. The full nonlinear 
band structure of the system is then determined for a 
two-color potential. We show that period-doubled states 
in the usual one-color lattice are directly connected with 
Bloch waves in the two-color case. Period-doubled states 
have a periodicity which is twice that of the underlying 
lattice. This allows for a novel interpretation of swallow- 
tails, a key physical property of the band structure of a 
BEC. The study of period-doubled states has been impor- 
tant in other systems, for instance in optics [rTHT^ |. be- 
cause it offers experimental access to the period-doubling 
route to chaos. 

The band structure of a BEC in an optical lattice 
has previously been studied both analytically and nu- 
merically. An array of techniques have been devel- 
oped. The sinusoidal potential, which is the experi- 
mental case, has been investigated by expanding the 
wavefunction in a Fourier series and minimizing the en- 
ergy d Q E El E3- Two other potentials have 
been studied in detail. A Jacobi elliptic potential, which 
asymptotically approaches the sinusoidal case, has the 
advantage that it emits a useful class of exact solu- 
tions [Tsl . Eii . l20l ]. On the other hand, the Kronig-Penney 



potential can be solved completely and exactly, as we 
have shown in our previous work |2l|, l22l I23I ] . We com- 
pare these three models with the aid of the Bloch wave 
representation. For the Jacobi elliptic potential this re- 
quires a reinterpretation of the original result Ell- in 
addition, wc explain the longstanding open problem of 
why certain density offsets in this class of solutions lead 
to instability. Wc find that the form of the potential has 
no effect on the superfluid and other physical properties 
of the system. Of the closed-form analytical methods 
available for the three potentials, only the Kronig-Penney 
potential admits a full description of the band structure. 

The two-color lattice has already received some atten- 
tion in the literature. Roth and Burnett showed that 
when the period of the second lattice is much longer 
than the first, so that it forms an envelope, new quantum 
phases are introduced into the quantum problem |24| . We 
will consider the superfluid phase only, which is obtained 
when the lattice height is on the order of or smaller than 
the chemical potential. Very recently, Louis et ai. have 
studied gap solitons in a two-color lattice for which the 
second lattice has half the period of the first [2j|. We 
study the same potential, but for Bloch waves, rather 
than gap solitons. Bloch waves have the same period as 
the lattice, whereas gap solitons are envelope solutions. 
Using the Kronig-Penney potential, we link the period- 
doubled states of the one-color lattice to Bloch waves 
in the two-color case. Adiabatically turning off one of 
the frequencies in a two-color lattice, one can clearly ob- 
serve the origin of swallowtails. Unlike previous expla- 
nations of this intriguing property of nonlinear periodic 
systems |26j . ours encompasses both repulsive and at- 
tractive BEC's. 

Experiments on BEC's in optical lattice potentials pro- 
ceed as follows. Alkali-metal bosonic atoms are cooled to 
the quantum degenerate regime. The interference pat- 
tern of two counter propagating lasers is used to cre- 
ate a sinusoidal potential shift caused by the ac Stark 
effect induced by the dipole interaction with the laser 
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field on the atoms' center of mass motion [271 l2Sj . A 
small frequency detuning of one of the lasers allows for 
examination of the different quasimomentum states and 
their stability properties by creating a traveling wave in- 
terference pattern moving at the velocity v = (X/2)Sv, 
where A is the wave leng th of the first beam and 8v is 
the detuning j2£J, l30l l31j . The two-color lattice can then 
be created by the superposition of two lasers with fre- 
quencies which differ by a factor of two. This can be 
achieved using second harmonic generation with nonlin- 
ear crystals. Complex lattice configurations with several 
lasers have already been performed experimentally |32j |. 
We present a dynamical study of key observables for the 
two-color lattice, such as the instability time, that can be 
experimentally investigated with minor changes to cur- 
rent apparatus. There have been many theoretical stud- 
ies of stability properties of BEC's in one-color lattices 
(see [ll[ll|Ij|M0IMIIl,to name a few), though 
not of the two-color case. We take advantage of the ana- 
lytical simplicity of the Kronig-Penney potential in order 
to correctly seed our numerical studies of the mean field. 

Thus, the three main objectives of this article are the 
following. First, in Sees. ITT1 and IIIII the qualitative prop- 
erties of a condensate in a periodic potential are shown 
to be independent of the specific form of the potential. 
This justifies the use of a Kronig-Penney potential for the 
rest of our study. Second, in Sec. IIV1 the nonlinear band 
structure of the two-color lattice is described fully and 
analytically. This is used to interpret the appearance of 
swallowtails in the one-color lattice in a novel way. Third, 
in Scc.0 numerical simulations are used to make a con- 
nection between the properties of the two-color lattice 
and current experiments on BEC's in one-color lattices. 
Finally, concluding remarks are made in Sec. IVII 



II. 



THE NONLINEAR SCHRODINGER 
EQUATION AND BLOCH WAVES 



The stationary states of the mean field of a BEC in 
an external periodic potential V(x) are governed by the 
nonlinear Schrodingcr equation |38l l39| , also called 
the Gross-Pitaevskii equation |4C 
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(1) 



where p is the eigenvalue and g characterizes the effec- 
tive quasi-onc-dimcnsional two-body atomic interaction. 
A quasi-one dimensional condensate can be created by 
making an elongated trap that is tightly bound in the 
transverse directions. In Eq. Q, length has been scaled 
by the lattice spacing d and energy has been scaled by 
2E /tt 2 where 



En 



2md 2 



(2) 



is the recoil energy, or kinetic energy of a particle at the 
edge of the first Brillouin zone, and m is the particle's 



mass. With this scaling, the renormalized ID coupling 
is g = 2a s u>md/h, where harmonic confinement in the 
transverse directions has been assumed with frequency 
u), and a s is the two body s-wave scattering length. 
The wavefunction ^ is given by 



$ = y/ p(x) exp[i(j)(x) — ifit] 



(3) 



where p(x) is the linear density of the condensate and 
the local superfluid velocity is given by v = d(f>/dx. Sub- 
stituting Eq. m t° Eq. QJ, one finds that the phase is 
related to the density by 



cj>{x) 



x dx' 

o pO') 



(4) 



where a is a real constant of integration. Note that we 
have chosen </>(0) = 0. 

When describing the band structure of a system, one 
constructs the solutions in the usual Bloch form, 



(5) 



where q is the quasimomentum and f q (x) = f q (x+l), i.e., 
has the same period as the lattice. The band structure 
E(q) can be determined from the relationship between 
the quasimomentum and the energy per particle in terms 
of the density and phase of the solution: 



1 dx' 



(6) 



- = - / 1 d a; '(i|a aif *| 2 + f|*|* + v( a /)|*| a ), (?) 

n n Jq 2 2, 

where E/n is the energy per atom, q is the quasimomen- 
tum, and n is the number of atoms per lattice site. The 
latter is given by the normalization 



dx'p(x') 



(8) 



The quasimomentum is the phase jump between adja- 
cent sites of the lattice and is proportional to the mean 
superfluid velocity of the system. We scale the quasimo- 
mentum by 



q a = w/d, 



(9) 



which is the quasimomentum at the edge of the Brillouin 
zone. 

The band structure is modified from the well-known 
solutions of the linear Schrodingcr equation with a peri- 
odic potential. There are two physical regimes [2lJ. For 
\g\n < Vb, which we term the regime of weak nonlinearity, 
the linear band structure is simply perturbed up or down, 
depending on the sign of g. For \g\n 3> Vb, which cor- 
responds to the regime of strong nonlinearity, the bands 
wrap back around on themselves to form loop structures, 
or swallowtails. These swallowtails have previously been 
described in terms of the superfluid screening properties 
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of the condensate [26j. This explanation is only valid 
for repulsive condensates, despite the fact that swallow- 
tails also appear in the attractive case. This motivates 
the need for another explanation, which we provide in 
Sec. El 

There have been many methods used to determine the 
nonlinear band structure of BEC's in periodic structures. 
In the following section, three analytical methods arc de- 
scribed and evaluated. In addition to these analytical ap- 
proaches, many numerical methods have been developed. 
These include Hamiltonian perturbation theory |42| and 
the study of accelerating lattices 01 ■ We consider only 
analytical methods of obtaining stationary states. 



III. ANALYTICAL METHODS IN NONLINEAR 
BAND THEORY 

Of the many methods used in nonlinear band theory, 
analytical methods allow for greater flexibility in describ- 
ing solution types, as they result in general expressions 
which describe all parameter regimes simultaneously. In 
this section, wc compare three methods that use differ- 
ent periodic potentials and interpret their solutions in the 
Bloch wave representation. It is found that they all pro- 
duce the same band structure for the lowest band. Hence, 
the exact form of the potential is unimportant. Never- 
theless, only the method based on the Kronig-Penney 
potential is able to describe all bands analytically since 
an exact solution method exists for piecewise constant 
potentials [USQ- 



A. Solution by Cancellation 

The first analytical method for obtaining solutions to 
the NLS with a periodic potential in the context of a BEC 
was given by Bronski, Carr, Deconinck, and Kutz for 
a Jacobi elliptic potential of form 



class of solutions with a period equal to that of the lat- 
tice, where the density was assumed to be of the form 



V(x) = -V {1 - 2sn 2 [2K(k)x, k}} , 



(10) 



where sn is one of the Jacobi elliptic functions [45| • The 
Jacobi elliptic functions are generalized periodic func- 
tions characterized by an elliptic parameter k € [0,1]. 
They approach circular and hyperbolic trigonometric 
functions as k — * and k — » 1, respectively. For k not 
exponentially close to unity, the potential is similar to 
the sinusoidal case. As k approaches unity, the period of 
the lattice becomes much greater than the width of the 
localized variations of the potential, thereby approaching 
a lattice of delta functions, i.e., a Kronig-Penney poten- 
tial. The period of the lattice is 2K(k) £ [tt, oo] where 
K(k) is a complete elliptic integral of the first kind [i^ . 

Bronski et al. were able to show that it is possible 
to choose a suitable ansatz for the density such that the 
nonlinear term cancels the potential term, thereby leav- 
ing a system that is effectively free. They found a special 



p(x) = Asn 2 (bx, k) + B . 
The following conditions must then be met: 
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a 



b 2 k 2 + 2V 
9 

B(B + A)(b 2 +Bg- 2 -^). 



1 



/i = -(b z + Ag + 3Bg)- 



A 
BVp 
A 



(11) 

(12) 

(13) 
(14) 



Machholm et al. [16j showed that this solution corre- 
sponds to the edge of the Brillouin zone for the exactly 
sinusoidal case, k = 0. If the parameters of the potential, 
Vb, fc, and 6, and the interaction strength gn are set, the 
wavefunction is completely determined without any free 
parameters. Nevertheless, the solutions are of more gen- 
eral use. In fact, although not described in this way by 
the original authors, this family of solutions can be used 
to map out the whole lowest energy band. It is possible 
to change the elliptic parameter k to get a spectrum of 
solutions. Since the Jacobi elliptic functions closely ap- 
proximate the trigonometric functions for all k except ex- 
ponentially close to unity, there is a wide range of values 
for k where the potential is approximately sinusoidal. In 
addition, the sign of the potential coefficient Vb does not 
significantly change the form of the potential and so this 
is another parameter that can be changed to determine 
the energy bands. Therefore, when the elliptic parameter 
is varied, and the solutions determined via Eqs. (|TT1) - ((T4*jl . 
the complete lowest energy band can be analytically de- 
termined for arbitrary interaction strength. 

An example of the band structure for strong nonlin- 
carity is shown in Fig. ^a). The interaction strength 
is chosen to be a factor of ten larger than the strength 
of the potential, causing the appearance of swallowtails. 
Only the lowest energy band can be extracted with this 
class of solutions. This shows that these simple exact 
solutions are sufficient to describe measurable properties 
of the condensate such as breakdown of superfluidity for 
critical values of the nonlineari tv 129.1311 . 

In the work of Bronski et al. [18L fl9l | , the linear stabil- 
ity properties of the solutions were proved for the case of 
constant, or trivial phase. No such proof was discovered 
for the case of non-trivial phase solutions. The latter are, 
for example, the only stable ones in the exactly sinusoidal 
case of k = 0. It was found numerically that non-trivial 
phase solutions for repulsive condensates with large off- 
sets B were stable, while for smaller B they were unsta- 
ble. No explanation for this stability property was found. 
After we recast the solutions in Bloch form and used them 
to determine the band structure, the stability properties 
become immediately apparent. Solutions on the upper 
edge of the swallowtail are known to be unstable, as they 
represent an energy maximum |l6l l2ll l26j ; these corre- 
spond to small B. Solutions on the main part of the 
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FIG. 1: Band structure for a strongly repulsive condensate, gn — 10, in a periodic potential. Shown are results obtained via 
the analytical methods of (a) Bronski et al. [Tsll with a Jacobi elliptic potential (\Vo\ = 1), (b) Machholm et al. flrjl with a 
sinusoidal potential (Vb = 4), and (c) Seaman et al. |2ll| with a Kronig- Penney potential (Vo = 1). 



band and the lower edge of the swallowtail represent an 
energy minimum and are stable; these correspond to large 
B. For instance, in Fig.QIa), the main part of the band 
from q = to q = 1 corresponds to B £ [0.559, 0.509]; the 
lower edge of the swallowtail near the edge of the Bril- 
louin zone corresponds to B £ [0.509, 0.277]; and the up- 
per edge, which is unstable, corresponds to B £ [0.277, 0]. 
We note that, in comparing with Bronski et al., they nor- 
malized the wavefunction to gn rather than having gn 
as the coefficient of the nonlinearity in the NLS. In their 
numerical studies of stability, they did not hold their nor- 
malization (our nonlinear coefficient) fixed. In Fig.^a), 
we fix the nonlinearity to be gn = 10|Vb|. 



B. Solution by Three Mode Approximation 

In a subsequent work by Machholm, Pethick and 
Smith |Tfij . a more complete method to determine the 
nonlinear band structure was introduced that used a si- 
nusoidal potential 

V(x) = 2Vbcos 2 (7re) . (15) 

The wavefunction was expanded in a Fourier series for 
a particular quasimomcntum and the energy was mini- 
mized. This method uses the exact physical form of the 
typical experimental lattice potential that a BEC is held 
in but requires the root finding of many free parameters. 
Machholm et al. showed that one can obtain analytical 
solutions by using a three-mode Fourier spectrum: 

= Vn~e lqx (a + aie i27ra + a^e" 12 ™) , (16) 



where ao, ai and a_i are real coefficients. Due to the 
normalization condition on the wavefunction, there are 
two real free parameters of the solution. 

Using the three-mode expansion of the wavefunction, 
it is possible to extract the first two bands and the lower 
part of the third band, as shown in Fig^b). Notice that, 
unlike with the previous method, full information from 
the first two bands can be extracted, instead of only the 
lowest band. However, the three-mode approximation 
overestimates the width of the swallowtails [ly . If more 
Fourier components arc included, for a total between five 
and ten, the first three bands can be described to within 
1% accuracy. For higher bands even more Fourier com- 
ponents are needed. This method was later extended to 
include period-doubled states [l7j . 

In addition to determining the band structure, Mach- 
holm et al. were able to determine several analytical 
expressions associated with the stability properties of 
the condensate. In particular, the interaction strength 
at which swallowtails first appear in the lowest band, 
and the width of the swallowtails for a given interaction 
strength were determined analytically. Studies of the sta- 
bility of the condensate to small perturbations were also 
performed. 

C. Solution via a Piece-wise Constant Lattice 

In comparison, in our previous work [2l| a lattice of 
delta functions, a Kronig-Penney potential, was used 

V(x)=V S(x-j). (17) 



We presented the complete set of Bloch waves solutions in 
analytical form by solving the piece- wise-constant poten- 
tial case and using the appropriate boundary conditions 
to fix the parameters of the solution. The potential, how- 
ever, is quite different from the experimentally created 
sinusoidal lattice. The optical lattice potential is com- 
posed of a single Fourier component while the Kronig- 
Penney potential is a comb of equally weighted Fourier 
components. We proved that the most general form of 
the density over a finite interval of constant potential is 
given by 

k 2 b 2 

p(x) = B H sn 2 (6x + x ,k) , (18) 

9 

where again sn is one of the Jacobi elliptic functions. It 
should be noted that the square of any Jacobi elliptic 
function can be related linearly with the square of ev- 
ery other Jacobi elliptic function. Thus all twelve Jacobi 
elliptic functions are possible solutions. The chemical po- 
tential \i and phase prefactor a (see Eq. |@J) are given 
by 

fx=^[b 2 (l + k 2 )+3Bg], (19) 
a 2 = B(k 2 b 2 /g + B){b 2 + Bg). (20) 

The energy bands can be determined by varying one of 
the parameters in the wavefunction, such as B or b, and 
determining the quasimomcntum and the energy from 
Eqs. I© and J7J. In particular, the parameter scaling b 
was varied since it is closely related to the quasimomen- 
tum. The offset B and elliptic parameter k are deter- 
mined by number conservation and the boundary condi- 
tions across the delta functions. 

The first three bands for the Kronig-Penney lattice 
are presented in Fig.^c). Although only the first three 
bands are presented, our method can be used to deter- 
mine the higher bands with no additional computational 
intensity, unlike with the technique of Sec. IIII Bl For all 
energy bands, all that is required computationally is the 
root finding of two parameters, the offset B and the ellip- 
tic parameter k. The stability properties of the conden- 
sate were then numerically determined by dynamically 
evolving an initial state composed of the stationary so- 
lution plus a small amount of white noise. It was found 
that the stability of the energy bands depends on the 
interaction strength. For weak nonlinearity, the bands 
become unstable for quasimomcntum greater than ap- 
proximately 7r/2. For strong nonlinearity, the tops of the 
swallowtails are unstable but the remainder of the first 
band is stable. The second and higher bands are always 
unstable. 

In spite of the differences in the potentials, all three 
methods find a similar nonlinear band structure. In ad- 
dition, the stability properties of the three models were 
also the same. Therefore, to study the main physical as- 
pects of this system, any method can be used. All three 
methods are useful in determining the first band. For the 
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FIG. 2: A sketch of the potential is illustrated for the three 
cases of Figs. [3] and [I] with appropriately normalized boxes 
representing delta functions: (a) a one-color lattice of period 
d and A = 0; (b) a two-color lattice with a small difference 
in the potential strengths, period 2d and A = 0.2; (c) a one- 
color lattice with period 2d and A = 1. Adiabatically tuning 
the system from (a) to (c) allows for a novel explanation of 
the appearance of swallowtails in the band structure. 



full band structure, only our method is analytically and 
numerically tractable. 

IV. THE TWO-COLOR LATTICE AND 
FORMATION OF SWALLOWTAILS 

Current experiments with BEC's utilize only one-color 
lattices in the weakly nonlinear regime. Therefore there 
has been no experimental observation of the swallowtail 
structure, which requires strong nonlinearity. Two-color 
lattices, on the other hand, allow one to produce swallow- 
tails even for weak nonlinearity, as we shall show in the 
following. Moreover, adiabatically tuning the two-color 
lattice allows for an explanation of swallowtails for both 
repulsive and attractive BEC's. 

A two-color lattice is produced by adding a second fre- 
quency component of twice the fundamental frequency, 
or half the period. Such a straightforward modification 
of existing experiments can be made via second har- 
monic generation with nonlinear crystals. In Sec. 11111 we 
showed that the exact form of the potential is unimpor- 
tant. Therefore, to study the two-color lattice, we use a 
Kronig-Penney-like potential with delta functions of two 
different strengths: 

V(x) = V J2 [(l-A)<y(x-2j) + (l + A)*(x-2j+l)], 

j — — oo 

(21) 

where A £ [0,1]. The strengths of the two sublatticcs 
are (1 — A)Vq and (1 + A)Vb, respectively. The form 
of Eq. (|21|l allows the potential to be tuned continu- 
ously from the one-color lattice of strength Vo and pe- 
riod d (A = 0) to the one-color lattice of strength 2Vb 
and period 2d (A = 1). This holds the average poten- 
tial strength constant. For < A < 1 the potential is 
two-color. Figure |3 sketches the progression between the 



two lattices, where the delta functions are represented by 
square functions of fixed width and varying heights. We 
will show how period-doubled solutions of the A = lat- 
tice map onto Bloch-wave solutions of the A = 1 lattice, 
including the swallowtails. 

The stationary states of the NLS with the potential of 
Eq. H21J) are determined in the same manner as was de- 
scribed in Sec. lIII 01 with the addition of a second bound- 
ary condition. Period-doubled solutions can be obtained 
via the transformation ^(x) = f q (x) exp(i2qx), while 
Bloch solutions were found via ty(x) = f q (x) exp(iqx) . 
In general, one can obtain solutions of the form ty(x) = 
fq{x) exp(imgx) , with m an integer; we focus on the 
period-doubled case, m = 2, and the Bloch wave case, 
m = 1. 

We first illustrate period-doubled solutions for the lat- 
tice of Fig. OJa), A = 0. There are two types of period- 
doubled solutions, trivial and non-trivial. The triv- 
ial period-doubled solutions simply reproduce the Bloch 
waves. That is to say, solutions of period d are also 
trivially solutions of period 2d. These arc illustrated in 
Fig. Ela) as thin curves for the case of weak nonlincar- 
ity Period-doubled solutions extend across a Brillouin 
zone of half the quasimomcntum-domain as that of Bloch 
waves. Thus the domain of Fig. El (? G [0,1]) is half 
that of Fig. Q] (<7 <= [ — F !])■ The bands are required to 
be symmetric around the center of the Brillouin zone. 
One reflects the band in the right half of the Bloch Bril- 
louin zone around its center at q = 0.5 to obtain the 
trivial period-doubled solutions. This leads to the two 
branches which make the form of an 'x' in Fig. OJa) (thin 
black curves). The non-trivial period-doubled solutions 
are shown as a thick blue curve in Fig.|2{a). These form a 
saddle between the trivial solution branches. The form of 
the density for the non-trivial period-doubled solutions is 
shown as the solid curve in Fig. [5] The non-trivial solu- 
tions are caused by the nonlinearity in Eq. Q ; they have 
no analog in the linear Schrodinger equation. Comparing 
Fig. |2| to the swallowtails illustrated in Fig. ^ we make 
the following conjecture: non-trivial period-doubled so- 
lutions for a lattice of period d appear as the saddle of 
the swallowtail for Bloch waves for a lattice of period 2d. 

This conjecture is supported by tuning A from zero 
to unity. We illustrate this tuning in the subsequent 
panels of Fig. Figure |^b) shows the case A < 1. 
Figure |3fc) shows the endpoint, A = 1. The shape of 
the lattice in the panels (a)-(c) is sketched in Fig. [21 In 
Fig. Otb) , the trivial and non-trivial period-doubled so- 
lutions of Fig. |3{a) separate into two bands. The lower 
band (solid curve) consists of the lower part of the triv- 
ial period-doubled solutions, with the non-trivial period- 
doubled solutions forming the upper edge of the swal- 
lowtail. The upper band (dashed curve) consists of the 
upper part of the trivial period-doubled solutions. One 
observes that there is a small gap between the two bands. 
Moreover, unlike in the case of the one-color lattice, one 
obtains a swallowtail even for weak nonlinearity As A 
is increased, this gap increases. Figure |^c) shows the 
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FIG. 3: (color online) Band structure for a two-color lat- 
tice with weak nonlinearity, gn = Vo. (a) Shown are period- 
doubled solutions for a one-color lattice of period d and A = 
(Thin black curve: trivial solutions; thick blue curve: non- 
trivial); (b),(c) Bloch-wave solutions with period 2d (Solid 
curve: lowest band; dashed curve: second band), (b) A two- 
color lattice with A <C 1. Note that the central swallowtail is 
derived from the non-trivial period-doubled solution of panel 
(a), (c) A one-color lattice with A = 1. The swallowtail 
has disappeared due to the weak nonlinearity. The lattices 
associated with (a)-(c) are sketched in Fig. [21 
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FIG. 4: (color online) Same as Fig. but for strong nonlin- 
earity, gn = WVq. In (b), (c) note that the upper edge of the 
central swallowtail in the lower band is derived from the non- 
trivial period-doubled solution of panel (a) (blue thick curve). 
The lattices associated with (a)-(c) are sketched in Fig. [5] 



case A = 1 . The system is now again in a weakly nonlin- 
ear regime with a one-color lattice. Thus the swallowtail 
disappears. 

One observes that, in general, the first two Bloch-wave 
bands of a one-color lattice of period 2d (panel (c)) can 
be derived from the trivial period-doubled solutions of a 
lattice of period d (panel (a)) via a two-color intermediate 




q (units of q Q ) 

FIG. 5: The condensate density profiles are illustrated for 
three cases: A = 0, non-trivial period-doubled solutions 
(solid curve); A = 1/2, upper edge of the swallowtail (dashed 
curve); and A = 1, the lowest band (dotted curve). The pa- 
rameters are q = 0.5 and gn — Vo, which is the regime of 
weak nonlinearity, as in Fig. [3] 



lattice (panel (b)). This follows from the fact that the 
second Bloch-wavc band has two density peaks per site. 

In the strongly nonlinear regime, the swallowtail forms 
in a similar manner but does not disappear for A = 1. 
We illustrate this sequence in Fig. again for A = 0, 
A <C 1, and A = 1. One clearly sees that the upper edge 
of the swallowtail in the one-color lattice of period 2d 
is derived from the non-trivial period-doubled solutions 
of the one-color lattice of period d. The lower edge of 
the swallowtail is derived from the trivial period-doubled 
solutions. 

This introduces a novel perspective on the origin of the 
swallowtails in the nonlinear band structure. In the fig- 
ures we illustrated the case of a repulsive condensate, as 
that is the most common experimentally. The same ar- 
gument holds for attractive condensates, where the swal- 
lowtails form below, rather than above the bands. In 
particular, in our study of the one-color Kronig-Pcnncy 
potential |2I|. we found it intriguing that no swallowtail 
formed on the lowest band for the attractive case. 

Previous interpretations of swallowtails were based on 
the superfluid screening properties of the condensate |26j | . 
In this argument, the condensate sees the quadratic free 
particle dispersion up to the sound speed. When written 
in the form of the Bloch ansatz, these quadratic curves 
repeat in each Brillouin zone. When the sound speed is 
such that the maximum quasimomcntum of each curve 
overlaps with the curve from the adjoining Brillouin zone, 
one obtains swallowtails. This does not apply to attrac- 
tive condensates, since the swallowtails form on the lower 
edge of the bands and not at all on the lowest band. 

In contrast, our argument based on period-doubling 
applies to both repulsive and attractive condensates. For 
attractive nonlinearity, the non-trivial period-doubled 



states form a convex saddle on the lower edge of the 
'x' of Fig. Ola), rather than the upper edge. As A is 
tuned from zero to unity, the upper part of the 'x' sepa- 
rates from the lower part, carrying the swallowtail with 
it. Thus no swallowtail can form on the lowest band. 

Thus, we understand swallowtails to be adiabatically 
connected to period-doubled solutions of a lattice of half 
the period. They originate in non-trivial period-doubling 
brought about by the nonlinear term in the NLS. 

Finally, since the density is a key observable in experi- 
ments on BEC's, in Fig. [5] density profiles are illustrated 
for weak nonlinearity and A = 0, 1/2, 1. The plot is 
made for quasimomentum q = 0.5 and the energy per 
particle such that the solution lies at the center of the 
non-trivial period-doubled solution (thick blue curve) in 
Fig-EIa), at the top of the central swallowtails on the low- 
est band in Fig.^b), and on the lower band in Fig.[3{c). 



V. STABILITY PROPERTIES AND 
COMPARISON WITH EXPERIMENTS 



0.75 




0.4 0.6 
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FIG. 6: Shown is the inverse of the instability time 7 = T^ 1 
of stationary states in the lowest band, as a function of the 
quasimomentum. White noise was added to the initial state 
and the wavefunction was evolved numerically for the NLS 
with a one-color lattice and weak nonlinearity, gn = Vq. 



Several important experiments have been performed 
with superfluid BEC's in optical lattices. In particular, 
an experiment by Fallani et al. [2^, ^(| studied the loss 
rate of the BEC as the quasimomentum was varied. This 
loss rate is expected to be monotonically related to the 
instability time of the condensate. This observable can 
be modelled theoretically by adding white noise to an ini- 
tial stationary state on the lattice in simulations. In this 
way, we investigate the nonlinear stability properties of 
Bloch-wave solutions to the two-color lattice. Note that 
unstable solutions which have lifetimes longer than ex- 
perimental timcscales will appear experimentally stable. 

In our simulations, Eq. Q was dynamically evolved us- 
ing a variable step fourth-order Runga-Kutta algorithm 
in time and a pseudospectral method in space. We chose 
periodic boundary conditions in one dimension, with a 
sufficient number of sites so that the results were inde- 
pendent of the ring circumference. Two schemes for delta 
functions on a grid were considered: point defects, and 
narrow square potentials covering several grid points. For 
square potentials less than one tenth the lattice period 
d, differences between the two schemes were negligible. 
Simulations were performed over time scales on the order 
of those relevant to the experiments, namely, hundreds 
of milliseconds. The lattice spacing d is taken to be 1 fxm 
and the atomic mass to be that of 87 Rb. 

The stability of the system is quantified by the variance 
in the Fourier spectrum, 



a(t) 



£(/(p,*)-/(p,Q)) 2 

2£(/(P,0)) 2 



(22) 



where f(p,t) is the Fourier component of the wavefunc- 
tion at momentum p and time t. The variance approaches 
unity when the wavefunction has no component in com- 
mon with the initial state and zero when then the wave- 
function it is identical to the initial state. We define 



the instability time Tj to occur when the variance has 
increased to a = 0.5. 

A clear experimental observable is the rate of loss of 
atoms from the system, since the density of the sample 
can be imaged at successive times. This loss rate is ex- 
pected to be monotonically related to the inverse time in 
which instabilities arise. Therefore, we present two rel- 
evant studies of the instability times of the condensate. 
First, the growth rate of the instability, 



7 = 1/Ti , 



(23) 



is determined for a one-color lattice and weak nonlin- 
earity. This relates directly to the experiment of Fal- 
lani et at. Second, the lowest positive quasimomentum 
for which the system becomes unstable is evaluated as a 
function of A for the two-color lattice. 

Figure H3 presents the results of the first study. The 
solutions are stable until the quasimomentum where the 
nontrivial and trivial period-doubled solutions connect in 
the band structure (see Fig. Of a), as well as Ref. ^In- 
stability can also be determined from the effective mass 
of the system, 



1 



d 2 E/d 2 q 



(24) 



If the effective mass becomes negative, the NLS becomes 
effectively attractive. Since the ground state of an at- 
tractive condensate is localized, solutions with negative 
effective mass m* < are unstable. The point where the 
nontrivial and trivial period-doubled solutions connect 
corresponds to a change in the sign of the effective mass 
from positive to negative. Our figure qualitatively repli- 
cates the results found by Fallani et al. |3lJ (see figures 
therein |3lj h 

Figure \7\ presents the second study. In particular, the 
smallest quasimomentum for which the periodic system 
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FIG. 7: The lowest positive quasimomentum for which the 
system becomes unstable is determined as a function of A as 
the one-color lattice of period d is tuned via a two-color lattice 
to the one-color lattice of period 2d. The result of simulations 
of the NLS with weak nonlinearity, gn = Vo, are shown. 



first becomes unstable in a time less than 100 ms is plot- 
ted. The two end points at A = and A = 1 arc 
easily determined, since those systems are one-color lat- 
tices. As per our discussion in the previous paragraph, 
the Bloch wave solutions become unstable at the point 
where the non-trivial period-doubled solutions connect to 
the Bloch-wave band. The two one-color lattice systems 
have the same ratio of nonlinearity to potential strength, 
since for A = the potential strength is Vq and the 
number of atoms is n, while for A = 1 they are 2Vo and 
2n. Thus the existence of the non-trivial period-doubled 
states scales with the length of the Brillioun zone. The 
absolute quasimomentum where instability occurs then 
differs by a factor of two, since for a lattice of period 2c? 
the Brillouin zone is half the length as for d. In Fig. [7] this 
is q ~ 0.6 for A = and q ~ 0.3 for A = 1. In the case 
of the linear Schrodingcr equation, the endpoints would 
be 0.5 and 0.25, respectively. Note that our simulations 
are accurate to within a few percent. 

The intermediate points in Fig. [7| where < A < 1 
and one obtains a two-color lattice, can be understood 
by again considering where the effective mass becomes 
negative. As the energy band from the A = lattice sep- 
arates into the first two bands of the A = 1 lattice, the 
upper edge of the swallowtail shrinks. Recall that states 
on the lower edge of swallowtails are stable |2lJ ■ As A 
is increased, the swallowtail disappears, since the nonlin- 
earity is weak. The quasimomentum at which the band 
changes from concave up to concave down, i.e. where the 
effective mass changes sign, moves from right to left, as 
can be seen in the descending data points of Fig. [7| A 



similar argument holds for the case of strong nonlinearity. 
VI. CONCLUSION 

We have presented a novel way of understanding the 
appearance of swallowtails in nonlinear band structure. 
Swallowtails in a lattice of period 2c? are adiabatically 
connected to the non-trivial period-doubled solutions of 
a lattice of period d. By non-trivial we mean solutions 
not immediately obvious from the fact that a solution 
with period d is also a solution of period 2d. Non-trivial 
period-doubled solutions are caused by the nonlinearity 
in the mean field of the condensate [IT]]- They do not 
appear in the linear Schrodingcr equation and therefore 
have no analog in linear band theory. Our way of un- 
derstanding swallowtails is valid for both repulsive and 
attractive condensates, unlike previous explanations |26| . 

We used the two-color lattice to adiabatically connect 
the one-color lattice of period d and 2d. We showed that 
swallowtails appear in the band structure of the two-color 
lattice even for weak nonlinearity, unlike in the one-color 
case. We then made explicit predictions for the onset of 
instability in the band structure for the one- and two- 
color lattices based on numerical simulations. Instability 
is experimentally observable as an increased loss rate due 
to heating of the condensate [3lj , and weak nonlinearity 
is the present regime of experimental investigation. 

Our studies utilized a two-color Kronig-Penney poten- 
tial. In order to justify the use of this model we compared 
analytic solutions to the one-color Kronig-Penney |2lj . 
Jacobi elliptic [l^, and sinusoidal potentials |l6j |. We 
showed that each of these models resulted in a similar 
band structure. Moreover, by putting the exact solu- 
tions to the Jacobi elliptic potential in Bloch wave form, 
we were able to resolve a longstanding question concern- 
in g de nsity offsets and stability in the work of Bronski et 

al is, m. 

We note that there are other solution types than the 
ones we have considered. For instance, there are gap 
solitons [25l |3^, 23 and time-periodic solutions ^3. |48l. 
These can have novel features in the two-color lattice [25j . 
We emphasize that we have treated only the superfiuid 
phase of cold bosons on a lattice [49l liSCj . 
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